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Abstract. An algorithm for the simulation of the evolution of slightly entangled 
quantum states has been recently proposed as a tool to study time-dependent 
phenomena in one-dimensional quantum systems. Its key feature is a time- 
evolving block-decimation (TEBD) procedure to identify and dynamically update 
the relevant, conveniently small subregion of the otherwise exponentially large 
Hilbert space. Potential applications of the TEBD algorithm are the simulation of 
time-dependent Hamiltonians, transport in quantum systems far from equilibrium 
and dissipative quantum mechanics. In this paper we translate the TEBD 
algorithm into the language of matrix product states in order to both highlight 
and exploit its resemblances to the widely used density-matrix renormalization- 
group (DMRG) algorithms. The TEBD algorithm, being based on updating 
a matrix product state in time, is very accessible to the DMRG community 
and it can be enhanced by using well-known DMRG techniques, for instance 
in the event of good quantum numbers. More importantly, we show how it 
can be simply incorporated into existing DMRG implementations to produce a 
remarkably effective and versatile "adaptive time-dependent DMRG" variant, that 
we also test and compare to previous proposals. 
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1. Introduction 

Over many decades the description of the physical properties of low-dimensional 
strongly correlated quantum systems has been one of the major tasks in theoretical 
condensed matter physics. Generically, this task is complicated by the strong quantum 
fluctuations present in such systems which are usually modelled by minimal-model 
Hubbard or Heisenberg-style Hamiltonians. Despite the apparent simplicity of these 
Hamiltonians, few analytically exact solutions are available and most analytical 
approximations remain uncontrolled. Hence, numerical approaches have always been 
of particular interest, among them exact diagonalization and quantum Monte Carlo. 

Decisive progress in the description of the low-energy equilibrium properties 
of one-dimensional strongly correlated quantum Hamiltonians was achieved by the 
invention of the density- matrix renormalization-group (DMRG) Ej ■ It is concerned 
with the iterative decimation of the Hilbert space of a growing quantum system such 
that some quantum state, say the ground state, is approximated in that restricted 
space with a maximum of overlap with the true state. Let the quantum state of a 
one-dimensional system be 

i j 

where we consider a partition of the system into two blocks S and E, and where {\i)} 
and are orthonormal bases of S and E respectively. Then the DMRG decimation 
procedure consists of projecting on the Hilbert spaces for S and E spanned by 
the M eigenvectors and corresponding to the largest eigenvalues of the 
reduced density matrices 

/i5-TrB|^)(^| p£ = Tr5|V^)(V|, (2) 

such that Ps\Wa) = ^^W^a) ^'^d PE\Wa) = XV\^a)- That both density matrices have 
the same eigenvalue spectrum is reflected in the guaranteed existence of the so-called 
Schmidt decomposition of the wave function 

|V^) = ^AJi«f)K), A„>0, (3) 

a 

where the number of positive A^ is bounded by the dimension of the smaller of the 
bases of S and E. 

Recently gl 13 El E], the ability of the DMRG decimation procedure to 
preserve the entanglement of \4>) between S and E has been studied in the context of 
quantum information science pillUj. This blooming field of research, bridging between 
quantum physics, computer science and information theory, offers a novel conceptual 
framework for the study of quantum many-body systems piBllSll^fTllFll^ llUllllLlT^ 
[l3. 14 15, 16, 17 . New insights into old quantum many-body problems can be gained 
from the perspective of quantum information science, mainly through its eagerness 
to characterize quantum correlations. As an example, a better understanding of the 
reasons of the breakdown of the DMRG in two-dimensional systems has been obtained 
in terms of the growth of bipartite entanglement in such systems |3E|. 

More specifically, in quantum information the entanglement of between S and 
E is quantified by the von Neumann entropy of ps (equivalently, of pe), 

5(/5s) = -E^'l°g2AL (4) 
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a quantity that imposes a useful (information theoretical) bound M > 2*^ on the 
minimal number M of states to be kept during the DMRG decimation process if 
the truncated state is to be similar to On the other hand, arguments from 

field theory imply that, at zero temperature, strongly correlated quantum systems 
are in some sense only slightly entangled in d = 1 dimensions but significantly more 
entangled in ci > 1 dimensions: In particular, in d = 1 a block corresponding to I 
sites of a gapped infinite-length chain has an entropy Si that stays finite even in the 
thermodynamical limit / — > oo, while at criticality Si only grows logarithmically with 
I. It is this saturation or, at most, moderate growth of Si that ultimately accounts for 
the succeess of DMRG ind — 1. Instead, in the general d-dimensional case the entropy 
of bipartite entanglement for a block of linear dimension / scales as Si ^ Thus, 
in d = 2 dimensions the DMRG algorithm should keep a number M of states that 
grows exponentially with I, and the simulation becomes inefficient for large I (while 
still feasible for small 

While DMRG has yielded an enormous wealth of information on the static and 
dynamic equilibrium properties of one-dimensional svstems|18l I19| and is arguably 
the most powerful method in the field, only few attempts have been made so far to 
determine the time evolution of the states of such systems, notably in a seminal paper 
by Cazalilla and Marston 20^- This question is of relevance in the context of the 
time-dependent Hamiltonians realized e.g. in cold atoms in optical lattices |21l I22j . 
in systems far from equilibrium in quantum transport, or in dissipative quantum 
mechanics. However, in another example of how quantum information science can 
contribute to the study of quantum many-body physics, one of us (G.V.) has recently 
developed an algorithm for the simulation of slightly entangled quantum computations 
| 23| that can be used to simulate time evolutions of one-dimensional systems jl7) . 

This new algorithm, henceforth referred to as the time-evolving block decimation 
(TEBD) algorithm, considers a small, dynamically updated subspace of the blocks S 
and E in Eq. to efficiently represent the state of the system, as we will review 
in detail below. It was originally developed in order to show that a large amount 
of entanglement is necessary in quantum computations, the rationale there being 
quite simple: any quantum evolution (e.g. a quantum computation) involving only 
a "sufficiently restricted" amount of entanglement can be efficiently simulated in a 
classical computer using the TEBD algorithm; therefore, from an algorithmical point 
of view, any such quantum evolution is not more powerful than a classical computation. 

Regardless of the implications for computer science, the above connection between 
the amount of entanglement and the complexity of simulating quantum systems is of 
obvious practical interest in condensed matter physics since, for instance, in d — 1 
dimensions the entanglement of most quantum systems happens to be "sufficiently 
restricted" precisely in the sense required for the TEBD algorithm to yield an 
efficient simulation. In particular, the algorithm has already been implemented 
and tested successfully on spin chains J7], the Bose-Hubbard model and single-atom 
transistors and dissipative systems at finite temperature pS). 

A primary aim of this paper is to reexpress the TEBD algorithm in a language 
more familiar to the DMRG community than the one originally used in Refs. |17ll2;^j . 
which made substantial use of the quantum information parlance. This turns out to 
be a rewarding task since, as we show, the conceptual and formal similarities between 
the TEBD and DMRG are extensive. Both algorithms search for an approximation 
to the true wave function within a restricted class of wave functions, which can be 
identified as matrix product states , and had also been previously proposed under 
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the name of finitely-correlated states|27|- Arguably, the big advantage of the TEBD 
algorithm relies on its flexibility to flow in time through the submanifold of matrix 
product states. Instead of considering time evolutions within some restricted subspace 
according to a fixed, projected, effective Hamiltonian, the TEBD algorithm updates 
a matrix product state in time using the bare Hamiltonian directly. Thus, in a sense, 
it is the Schrodinger equation that decides, at each time step, which are the relevant 
eigenvectors for S and E in Eq. (PJ, as opposed to having to select them from some 
relatively small, pre-selected subspace. 

A second goal of this paper is to show how the two algorithms can be integrated. 
The TEBD algorithm can be improved by considering well-known DMRG techniques, 
such as the handling of good quantum numbers. But most importantly, we will 
describe how the TEBD simulation algorithm can be incorporated into preexisting, 
quite widely used DMRG implementations, the so-called finite-system algorithmic 
using White's prediction algorithm 28 . The net result is an extremely powerful 
"adaptive time-dependent DMRG" algorithm, that we test and compare against 
previous proposals. 

The outline of this paper is as follows: In Section El we discuss the problems 
currently encountered in applying DMRG to the calculation of explicitly time- 
dependent quantum states. Section|3reviews the common language of matrix product 
states. We then express both the TEBD simulation algorithm (Sec. 0J and DMRG 
(Sec. Isj in this language, revealing where both methods coincide, where they differ 
and how they can be combined. In Section we then formulate the modifications 
to introduce the TEBD algorithm into standard DMRG to obtain the adaptive time- 
dependent DMRG, and Section 13 discusses an example application, concerning the 
quantum phase transition between a superfluid and a Mott-insulating state in a Bose- 
Hubbard model. To conclude, we discuss in Section|Slthe potential of the new DMRG 
variant. 

2. Simulation of time-dependent quantum phenomena using DMRG 

The first attempt to simulate the time evolution of quantum states using DMRG is 
due to Cazalilla and Marston 20'. After applying a standard DMRG calculation 
using the Hamiltonian H{t = 0) to obtain the ground state of the system at t = 0, 
ji/jo), the time-dependent Schrodinger equation is numerically integrated forward in 
time, building an effective Hcff{t) — i?cff(0) -I- Vcs{t), where iJefi(O) is taken as the 
Hamiltonian approximating H{0) in the truncated Hilbert space generated by DMRG. 
Vos{t) as an approximation to V{t) is built using the representations of operators in 
the block bases obtained in the standard DMRG calculation of the t = state. V{t) 
contains the changes in the Hamiltonian with respect to the starting Hamiltonian: 
H{t) — + V{t). The (effective) time-dependent Schrodinger equation reads 



away. If the evolution of the ground state is looked for, the initial condition is obviously 
to take I "0(0)) = IV'o) obtained by the preliminary DMRG nm. Forward integration 
can be carried out by step-size adaptive methods such as the Runge-Kutta integration 
based on the infinitesimal time evolution operator 




m+st)) = {i-iH{t)6t)m)) 



(6) 



Adaptive time- dependent DMRG 



5 



where we drop the subscript denoting that we are deahng with effective Hamihonians 
only. The algorithm used was a fourth-order adaptive size Runge-Kutta algorithm 

m- 

Sources of errors in this approach are twofold, due to the approximations involved 
in numerically carrying out the time evolution, and to the fact that all operators live 
on a truncated Hilbert space. 

For the systems studied we have obtained a conceptually simple improvement 
concerning the time evolution by replacing the explicitly non-unitary time-evolution 
of Eq. lO by the unitary Crank-Nicholson time evolution 

m^st))^l^m^m))- (7) 

l + {H{t)5t/2 

To implement the Crank-Nicholson time evolution efficiently we have used a (non- 
Hermitian) biconjugate gradient method to calculate the denominator of Eq. ((JJ). In 
fact, this modification ensures higher precision of correlators, and the occurence of 
asymmetries with respect to reflection in the results decreased. 

It should be noted, however, that for the Crank-Nicholson approach only lowest- 
order expansions of the time evolution operator exp(— i/fJt) have been taken; we have 
not pursued feasible higher-order expansions. 

As a testbed for time-dependent DMRG methods we use throughout this paper 
the time-dependent Bose-Hubbard Hamiltonian, 

HBH{t) = -J bl,b^ + h\h^^, + ^ E "'('^^ " (8) 
1=1 i=l 

where the (repulsive) onsite interaction [/ > is taken to be time-dependent. This 
model exhibits for commensurate filling a Kosterlitz-Thouless-like quantum phase 
transition from a superfluid phase for u < Uc (with u — U / J) to a Mott-insulating 
phase for u > Uc- We have studied a Bose-Hubbard model with L = 8 and open 
boundary conditions, total particle number = 8, J = 1, and instantaneous switching 
from Ui = 2 in the superfluid phase to U2 — 40 in the Mott phase at i = 0. We consider 
the nearest-neighbor correlation, a robust numerical quantity, between sites 2 and 3. 
Up to 8 bosons per site (i.e. A^gito = 9 states per site) were allowed to avoid cut-off 
effects in the bosonic occupation number in all calculations in this Section. All times 
in this paper are measured in units of h/ J or 1/ J, setting h=l. Comparing Runge- 
Kutta and Crank-Nicholson (with time steps oi St ^ 5 x 10^^) we found the latter 
to be numerically preferable; all static time-dependent DMRG calculations have been 
carried out using the latter approach. 

However, Hilbert space truncation is at the origin of more severe approximations. 
The key assumption underlying the approach of Cazalilla and Marston is that the 
effective static Hilbert space created in the preliminary DMRG run is sufficiently 
large that \ipit)) can be well approximated within that Hilbert space for all times, 
such that 

e(i)-l-|(V'(t)|Vcxact(i))l (9) 
remains small as t grows. This, in general, will only be true for relatively short times. 
A variety of modifications that should extend the reach of the static Hilbert space 
in time can be imagined. They typically rest on the DMRG practice of "targeting" 
several states: to construct the reduced density matrix used to determine the relevant 
Hilbert space states, one may carry out a partial trace over a mixture of a small 
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number of states such that the truncated Hilbert space is constructed so that all of 
those states are optimally approximated in the DMRG sense: 

Ps = Tr£;|^)(7/^| ps^ Tr£;^aj|V'i)(V'i|- (10) 

i 

A simple choice uses the targeting of i/"|-0o)i for n less than 10 or so, 
approximating the short-time evolution, which we have found to substantially improve 
the quality of results for non-adiabatic switching of Hamiltonian parameters in time: 
convergence in M is faster and more consistent with the new DMRG method (see 
below). 

Similarly, we have found that for adiabatic changes of Hamiltonian parameters 
results improve if one targets the ground states of both the initial and final 
Hamiltonian. These approaches are conceptually very similar to targeting not only 
|i/'o)i but also OlV'o) and some -ff"0|V'o); n = 1, 2, 3, . . . in Lanczos vector dynamics 
DMRG |3(JI \61\ . or real and imaginary part of [H — uu — Eq + 177)^^01-00) in correction 
vector dynamics DMRG pTl I32| to calculate Green's functions 

^^ol^^TT ^— -OlV-o). (11) 

H — u — bjo + iri 

To illustrate the previous approaches, we show results for the parameters of the 
Bose-Hubbard model discussed above. Time evolution is calculated in the Crank- 
Nicholson approach using a stepwidth (5t = 5 • 10~^ in time units of h/ J targeting (i) 
just the superfluid ground state |i/;o) for Ui = 2 (Fig. P), (ii) in addition to (i) also 
the Mott-insulating ground state ji/'o) for U2 = 40 and H{t > 0)|V'o) (Fig- EJ, (iii) in 
addition to (i) and (ii) also Hit > 0)^1 V'o) and H{t > Ofl^po) (Fig. EJ. 

Wc have used up to M = 200 states to obtain converged results (meaning that 
we could observe no difference between the results for M = 100 and M = 200) for 
t < i, corresponding to roughly 25 oscillations. The results for the cases (ii) and (iii) 
are almost converged for AI = 50, whereas (i) shows still crude deviations. 

A remarkable observation can be made if one compares the three M — 200 curves 
(Fig. 0J, which by standard DMRG procedure (and for lack of a better criterion) 
would be considered the final, converged outcome, both amongst each other or to the 
result of the new adaptive time-dependent DMRG algorithm which we are going to 
discuss below: result (i) is clearly not quantitatively correct beyond very short times, 
whereas result (ii) agrees very well with the new algorithm, and result (iii) agrees 
almost (beside some small deviations at i « 3) with result (ii) and the new algorithm. 
Therefore we see that for case (i) the criterion of convergence in M does not give 
a good control to determine if the obtained results are correct. This raises as well 
doubts about the reliability of this criterion for cases (ii) and (iii) . 

A more elaborate, but also much more time-consuming improvement still within 
the framework of a static Hilbert space was proposed by Luo, Xiang and Wang 33 , 3lj. 
Additional to the ground state they target a finite number of quantum states at various 
discrete times using a bootstrap procedure starting from the time evolution of smaller 
systems that are iteratively grown to the desired final size. 

The observation that even relatively robust numerical quantities such as nearest- 
neighbor correlations can be qualitatively and quantitatively improved by the 
additional targeting of states which merely share some fundamental characteristics 
with the true quantum state (as we will never reach the Mott-insulating ground state) 
or characterize only the very short-term time evolution indicates that it would be 
highly desirable to have a modified DMRG algorithm which, for each time t, selects 
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0.5 1 1.5 2 2.5 3 3.5 4 
time 



Figure 1. Time evolution of the real part of the nearest-neighbor correlations in 
a Bose-Hubard model with instantaneous change of interaction strength at t = 0: 
superfluid state targeting only. The different curves for different M are shifted. 



(ii) L=8, U1=2, U2=40, J=1 




0.5 1 1.5 2 2.5 3 3.5 4 
time 



Figure 2. Time evolution of the real part of the nearest-neighbor correlations in 
a Bose-Hubard model with instantaneous change of interaction strength at t = 0: 
targeting of the initial superfluid ground state, Mott insulating ground state and 
one time-evolution step. The different curves for different M are shifted. 
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(iii) L=8, U1=2, U2=40, J=1 




0.5 1 1.5 2 2.5 3 3.5 4 
time 



Figure 3. Time evolution of the real part of the nearest-neighbor correlations in 
a Bose-Hubard model with instantaneous change of interaction strength at t = 0: 
targeting of the initial superfluid ground state, Mott insulating ground state and 
three time-evolution steps. The different curves for different M are shifted. 

L=8, U1=2, U2=40, J=1 




0.5 1 1.5 2 2.5 3 3.5 4 
time 



Figure 4. Comparison of the three M = 200 Crank-Nicholson calculations to 
adaptive time-dependent DMRG at M = 50: we target (i) just the superfluid 
ground state Ii/jq) for Ui = 2 (Fig. 0, (ii) in addition to (i) also the Mott- 
insulating ground state Ii/^q) for U2 = 40 and H{t > 0)|i/'o) (Fig. EJ, (iii) in 
addition to (i) and (ii) also H(t > 0)2|i/)o) and H{t > 0)^|i/'o)- The different 
curves are shifted. 
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Hilbert spaces of dimension M such that \'(p{t)) is represented optimally in the DMRG 
sense, thus attaining at all times the typical DMRG precision for M retained states. 
The presentation of such an algorithm is the purpose of the following sections. 



3. Matrix product states 

As both the TEBD simulation algorithm and DMRG can be neatly expressed in the 
language of matrix product states, let us briefly review the properties of these states 
also known as finitely-correlated states (?7ll26j . 

Wc begin by considering a onc-dimensional system of size L, divided up into sites 
which each have a local Hilbert space, Tii. For simplicity we take the same dimension 
iVgito at all sites. In such a system a product state may be expressed as 

It) = |<7i)(g)|a2)®...® Ictl), (12) 
where \(Ji) denotes the local state on site i. We can express a general state of the 
whole system as 



Cl ,...,cr_L 

^^Vrrk). (13) 

(T 

This general state exists in the Hilbert space TL = Hi^i vjiih. dimension (A^site)'^- 
A matrix product state is now formed by only using a specific set of expansion 
coefficients ■i/'cr- Let us construct this set in the following. To do this we define 
operators j4j[(Ti] which correspond to a local basis state \ai) at site i of the original 
system, but which act on auxiliary spaces of dimension M, i.e., 

A;N=^Aj,^N|«)(/3|, (14) 

where \a) and \j3) are orthonormal basis states in auxiliary spaces. For visualization, 
we imagine the auxiliary state spaces to be located on the bonds next to site i. If we 
label the bond linking sites i and i + 1 by i, then we say that the states \j3) live on 
bond i and the states |a) on bond i ~ 1. The operators Ai[ai\ hence act as transfer 
operators past site i depending on the local state on site i. On the first and last site, 
which will need special attention later, this picture involves bonds and L to the left 
of site 1 and to the right of site L respectively. While these bonds have no physical 
meaning for open boundary conditions, they are identical and link sites 1 and L as 
one physical bond for periodic boundary conditions. There is no a priori significance 
to be attached to the states in the auxiliary state spaces. 

In general, operators corresponding to different sites can be different. If this is 
the case the resulting matrix product state to be introduced is referred to as a position 
dependent matrix product state. We also impose the condition 

]A,[a,]A\[a,]^l, (15) 



E- 



which we will see to be related to orthonormality properties of bases later. An 
unnormalized matrix product state in a form that will be found useful for Hamiltonians 
with open boundary conditions is now defined as 

l'^) = Ef<'^^in^^NI0H))k), (16) 

cr V i=i / 
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where \4>l) and are the left and right boundary states in the auxihary spaces 
on bonds and L. They act on the product of the operators Aj to produce scalar 
coefficients 

L 

^PcT ^ {MWM^.Ur) (17) 

i=l 

for the expansion of 

Several remarks are in order. It should be emphasized that the set of states 
obeying Eq. 1)16(1 is an (arbitrarily constructed) submanifold of the full boundary- 
condition independent Hilbert space of the quantum many-body problem on L 
sites that is hoped to yield good approximations to the true quantum states for 
Hamiltonians with open boundary conditions. If the dimension M of the auxiliary 
spaces is made sufficiently large then any general state of the system can, in 
principle, be represented exactly in this form (provided that and |0_r) are chosen 
appropriately), simply because the 0{NsiteLM'^) degrees of freedom to choose the 
expansion coefficients will exceed N^^^^. This is, of course, purely academic. The 
practical relevance of the matrix product states even for computationally manageable 
values of M is shown by the success of DMRG, which is known f551 155) to produce 
matrix product states of auxiliary state space dimension A/, in determining energies 
and correlators at very high precision for moderate values of M. In fact, some very 
important quantum states in one dimension, such as the valence-bond-solid (VBS) 
ground state of the Affleck-Kennedy-Lieb-Tasaki (AKLT) model [2313111211, can be 
described exactly by matrix product states using very small M (M = 2 for the AKLT 
model). 

Let us now formulate a Schmidt decomposition for matrix product states which 
can be done very easily. An unnormalized state \tp) of the matrix-product form of Eq. 
H16|l with auxiliary space dimension M can be written as 

M 
a=l 

where we have arbitrarily cut the chain into S on the left and E on the right with 

|<T^), (19) 



) = E 



bL\YlA,[a,]\a) 

ies 

and similarly {w^), where {la}} are the states spanning the auxiliary state space on 
the cut bond. Normalizing the states IV^), \w^) and \w^) we obtain the representation 
M 

|V>) = EAo|«^f)kf) (20) 

a=l 

where in Aq the factors resulting from the normalization are absorbed. The 
relationship to reduced density matrices is as detailed in Sec. ^ 



4. TEBD Simulation Algorithm 



Let us now express the TEBD simulation algorithm in the language of the previous 
section. In the original exposition of the algorithm (221 , one starts from a 
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Li/ • : • ^L-l I 

/ /+; 

Figure 5. Bipartitioning by cutting bond I between sites I and / + 1. 



representation of a quantum state where the coefficients for the states are decomposed 
as a product of tensors, 

- "S^ r[i]'^i \[i]r[2]T2 \[2]t.[3]<t3 . . .-rW'Tt /on 

Ql,...,Ql,_l 

It is of no immediate concern to us how the T and A tensors are constructed 
exphcitly for a given physical situation. Let us assume that they have been determined 
such that they approximate the true wave function close to the optimum obtainable 
within the class of wave functions having such coefficients; this is indeed possible as 
will be discussed below. There are, in fact, two ways of doing it, within the framework 
of DMRG (see below) , or by a continuous imaginary time evolution from some simple 
product state, as discussed in Ref. PT] . 

Let us once again attempt a visualization; the (diagonal) tensors X^^\ i — 
1, . . . ,L ~ 1 are associated with the bonds i, whereas F^, i = 2, . . . , L ~ 1 links 
(transfers) from bond i to bond i — 1 across site i. Note that at the boundaries 
{i = 1, L) the structure of the F is different, a point of importance in the following. 
The sums run over M states jaj) living in auxiliary state spaces on bond i. A priori, 
these states have no physical meaning here. 

The F and A tensors are constructed such that for an arbitrary cut of the system 
into a part Si of length I and a part -Bl-/ of length L — I ai bond I, the Schmidt 
decomposition for this bipartite splitting reads 

i^) = E^Si^S>Kro, (22) 



with 



and 



E E rW"^L'|---rS-Ul'^i)®---®h)' (23) 



ai,...,Qi_i cri.,...,(Ti 



Qi,...,Qi,_i ai^i,...,aL 

|a,+i)(8)---® Idi), (24) 

where \^) is normalized and the sets of and {\wa,''~'')} are orthonormal. This 
implies, for example, that 

E(^S)' = 1- (25) 

We can see that (leaving aside normalization considerations for the moment) 
this representation may be expressed as a matrix product state if we choose for 

^«N=Eo,/3^L/3NI«>(/3| 

Al,p[a,]=T%^^xf, (26) 
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except for i — 1, where we choose 

Ai^[ai] = /„rW^^AW, (27) 
and for i — L, where we choose 

^a/3[^L]=rm^-5/3. (28) 

The vectors fa and gp are normahsed vectors which must be chosen in conjunction 
with the boundary states and \4>r) so as to produce the expansion 121(1 from this 
choice of the Ai. SpecificaUy, we require 

|'/>l) = (29) 

a 

where \a) and \j3) are the states forming the same orthonormal basis in the auxihary 
spaces on bonds and L used to express A^^^. In typical implementations of the 
algorithm it is common to take fa — 9a ~ ^a,i- Throughout the rest of the article 
we take this as the definition for ga and /q, as this allows us to treat the operators 
on the boundary identically to the other operators for the purposes of the simulation 
protocol. For the same reason we define a vector Aq' = 5a,i- 

In the above expression we have grouped T and A such that the A reside on the 
right of the two bonds linked by F. There is another valid choice for the Ai, which will 
produce identical states in the original system, and essentially the same procedure for 
the algorithm. If we set 

^L/^N^ArHFM-, (31) 
except for i — 1, where we choose 

Alp[a,]^faTf^\ (32) 
and for i = L, where we choose 

^o;3K] - Al^-ilFW^-g^, (33) 

then the same choice of boundary states produces the correct coefficients. Here we 
have grouped F and A such that the A reside on the left of the two bonds linked by 
F. It is also important to note that any valid choice of fa and gp that produces the 
expansion (|21|) specifically excludes the use of periodic boundary conditions. While 
generalizations are feasible, they lead to a much more complicated formulation of the 
TEBD simulation algorithm and will not be pursued here. 

To conclude the identification of states, let us consider normalization issues. The 
condition (|15|l is indeed fulfilled for our choice of A^ffi], because we have from (|24|l 
for a splitting at / that 

= Y.<-.<^M\^i)®\<r)^ (34) 
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so that from the orthonormahty of the sets of states {\wa'' *' ^')}aLu {|o'i)}<T =i ^^'^ 

{k7^-')}^Li, 

<yi a/37 '^f- 

a/3 

= ^<5„^|a)(/3| =Z. (35) 

a/3 

Let us now consider the time evolution for a typical (possibly time-dependent) 
Hamiltonian in strongly correlated systems that contains only short-ranged 
interactions, for simplicity only nearest-neighbor interactions here: 

i odd j even 

Fi^i+i and Gj j+i are the local Hamiltonians on the odd bonds linking i and 1, and 
the even bonds linking j and j + \. While all F and G terms commute among each 
other, F and G terms do in general not commute if they share one site. Then the 
time evolution operator may be approximately represented by a (first order) Trotter 
expansion as 

^-iHSt ^ -Q g-iF...+i5t -Q g-iG,,, + i5t ^ 0{5t^), (37) 
i odd j even 

and the time evolution of the state can be computed by repeated application of the 
two-site time evolution operators exp(— iGjj-|-i(5i) and exjp[—iFi^i^i5t) . This is a well- 
known procedure in particular in Quantum Monte Carlo ^Uj where it serves to carry 
out imaginary time evolutions (checkerboard decomposition). 

The TEBD simulation algorithm now runs as follows 1 17| : 

(i) Perform the following two steps for all even bonds (order does not matter): 

(i) Apply exp(— iGi^j+ifSt) to \'4>{t)). For each local time update, a new wave 
function is obtained. The number of degrees of freedom on the "active" 
bond thereby increases, as will be detailed below. 

(ii) Carry out a Schmidt decomposition cutting this bond and retain as in 
DMRG only those M degrees of freedom with the highest weight in the 
decomposition. 

(ii) Repeat this two-step procedure for all odd bonds, applying exp(— ii^; ;_|_i(5t). 

(iii) This completes one Trotter time step. One may now evaluate expectation values 
at selected time steps, and continues the algorithm from step 1. 

Let us now consider the computational details, 
(i) Consider a local time evolution operator acting on bond Z, i.e. sites I and I -I- 1, for 
a state |'(/'). The Schmidt decomposition of j-^) after partitioning by cutting bond I 
reads 

M 

E^SiOK^-')- (38) 

ai = l 
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Using Eqs. (gSl), 123 and ijSlJl, wc find 

K:;)k^)km)kf,r/'^^')- (39) 

We note, that if we identify \wail\) and juia,^"^*'^^' ) with DMRG system and 
environment block states \'wf^^_^) and \'w^^^^^)^ we have a typical DMRG state for 
two blocks and two sites 

i'/^)-EEEE^' (40) 

with 

V'mi-i.iai + imi + i = E ^™T-UU-ia, hl^atmi + l I^'+l] • (^1) 

The local time evolution operator on site Z + 1 can be expanded as 

C^M+i- E E (42) 

and generates = where 

1^')= E E E 

This can also be written as 
where 

er^^,,, = A[^-\i 5: <_,„,h']<i,,hVjL^;-;;. (44) 

(ii) Now a new Schmidt decomposition identical to that in DMRG can be carried out 
for cutting once again bond /, there are now MNsHc states in each part of the 
system, leading to 

IV') = E aSI*5)K"-')- (45) 

In general the states and coefficients of the decomposition will have changed compared 
to the decomposition (|38|l previous to the time evolution, and hence they are adaptive. 
We indicate this by introducing a tilde for these states and coefficients. As in DMRG, 
if there are more than M non-zero eigenvalues, we now choose the M eigenvectors 
corresponding to the largest Aq, to use in these expressions. The error in the final 
state produced as a result is proportional to the sum of the magnitudes of the discarded 
eigenvalues. After normalization, to allow for the discarded weight, the state reads 

M 

1^') = E^Sl^^DK'"')- (46) 

ai=l 
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Figure 6. Typical two-block two-site setup of DMRG as used here. 



Note again that the states and coefficients in this superposition are in general different 
from those in Eq. (|38|l : we have now dropped the tildes again, as this superposition 
will be the starting point for the next time evolution (state adaption) step. As is done 
in DMRG, to obtain the Schmidt decomposition reduced density matrices are formed, 
e.g. 

P£-Trs|^')(^'l 

XI ki + l)ka, + l)(^"a; + J('^/+ll 

X ( E 0S"iU+i(0r_r.Vi) I • (4^) 

If we now diagonahse pE, we can read off the new values of ^^^^^^ [^i+i] because 
the eigenvectors obey 

E <^,,,h+i]k/+i)kf,:i""^')- (48) 

We also obtain the eigenvalues, (Aq,)^. Due to the asymmetric grouping of T and A 
into A discussed above, a short calculation shows that the new values for A^ai^iai Wi\ 
can be read off from the slightly more complicated expression 

The states |wf') are the normalized eigenvectors of ps formed in analogy to pE- 

The key point about the TEBD simulation algorithm is that a DMRG-style 
truncation to keep the most relevant density matrix eigenstates (or the maximum 
amount of entanglement) is carried out at each time step. This is in contrast with 
time-dependent DMRG methods up to now, where the basis states were chosen before 
the time evolution, and did not "adapt" to optimally represent the final state. 



5. DMRG and matrix-product states 

Typical normalized DMRG states for the combination of two blocks S and E and two 
single sites (Fig. O have the form 

IV-) = E EE E V'm,_i^,<T, + im, + i|w^,_i)|o-i)|crz+l)|w^,^J (50) 
rrii-i CTi crj + i mj + i 

which can be Schmidt decomposed as 

iV') = E^"liOi<)- (51) 

It has been known for a long timel 35[l5^ that a DMRG calculation retaining M 
block states produces M x M matrix-product states for \tp). Consider the reduced 
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basis transformation to obtain the states of DMRG block S that terminates on bond 
I from those of the block terminating on bond I — 1 and those on a single site I, 

{wi,m\wi,)=A'_,^,['Ji]^ (52) 



such that 



Ka,_,mMl]\wi,^,)®Wl)- (53) 



The reduced basis transformation matrices A;[(T;] automatically obey Eq. (|15|l . which 
here ensures that {Iw^j^}} is an orthonormal set provided {\w^-^^__^)} is one, too. We 
may now use Eq. (|53|l for a backward recursion to express \w^^__^) via \w^^_^) and so 
forth. There is a complication as the number of block states for very short blocks is 
less than M . For simplicity, we assume that AI is chosen such that we have exactly 
iVg^j, = M . If we stop the recursion at the shortest block of size N that has M states 
we obtain 

K)= E E (54) 

r?j.jy^j...mi_i ai-.-ai 

where we have boundary-site states on the first N sites indexed by = {ui . . . cr^}. 
Similarly, for the DMRG block E we have 

{w'ri,^^^m+i\wf^,) = Aj;[^^^Jcri+i], (55) 

such that (again having N boundary sites) a recursion gives 

i<)= E E 

• ..A^-f__^^^^^^[a^_^]\ai+^ . ..aL), (56) 

with boundary-site states on the last N sites indexed by nij^_fj = {(T^_jy_|_]^ . . . ctl}. 

A comparison with Eqs. (|16|l . H18|l and (|19|l shows that DMRG generates position- 
dependent M X M matrix-product states as block states for a reduced Hilbert space 
of M states; the auxiliary state space to a bond is given by the Hilbert space of the 
block at whose end the bond sits. This physical meaning attached to the auxiliary 
state spaces and the fact that for the shortest block the states can be labeled by good 
quantum numbers (if available) ensures through (|52|l and (|55|l that they carry good 
quantum numbers for all block sizes. The big advantage is that using good quantum 
numbers allows us to exclude a large amount of wave function coefficients as being 
0, drastically speeding up all calculations by at least one, and often two orders of 
magnitude. Moreover, as is well known, DMRG can be easily adapted to periodic 
boundary conditions, which is in principle also possible for the TEBD algorithm but 
cumbersome to implement. Fermionic degrees of freedom also present no specific 
problem, and in particular, there exists no negative sign problem of the kind that is 
present in Quantum Monte Carlo methods. 

The effect of the finite-system DMRG algorithmT is now to shift the two free 
sites through the chain, growing and shrinking the blocks S and E as illustrated in Fig. 
[71 At each step, the ground state is redetermined and a new Schmidt decomposition 
carried out in which the system is cut between the two free sites, leading to a new 
truncation and new reduced basis transformations (2 matrices A adjacent to this 
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Figure 7. Finite-system DMRG algorithm. Block growth and shrinkage. For 
the adaptive time-dependent DMRG, replace ground state optimization by local 
time evolution. 



bond). It is thus a sequence of local optimization steps of the wave function oriented 
towards an optimal representation of the ground state. Typically, after some "sweeps" 
of the free sites from left to right and back, physical quantities evaluated for this 
state converge. While comparison of DMRG results to exact results shows that one 
often comes extremely close to an optimal representation within the matrix state 
space (which justifies the usage of the DMRG algorithm to obtain them), it has been 
pointed out and numerically demonstrated HT] that finite-system DMRG results 
can be further improved and better matrix product states be produced by switching, 
after convergence is reached, from the S»»E scheme (with two free sites) to an S»E 
scheme and to carry out some more sweeps. This point is not pursued further here, 
it just serves to illustrate that finite-system DMRG for all practical purposes comes 
close to an optimal matrix product state, while not strictly reaching the optimum. 

As the actual decomposition and truncation procedure in DMRG and the TEBD 
simulation algorithm are identical, our proposal is to use the finite-system algorithm 
to carry out the sequence of local time evolutions (instead of, or after, optimizing 
the ground state), thus constructing by Schmidt decomposition and truncation new 
block states best adapted to a state at any given point in the time evolution (hence 
adaptive block states) as in the TEBD algorithm, while maintaining the computational 
efficiency of DMRG. To do this, one needs not only all reduced basis transformations, 
but also the wave function in a two-block two-site configuration such that the 
bond that is currently updated consists of the two free sites. This implies that \'4>) 
has to be transformed between different configurations. In finite-system DMRG such 
a transformation, which was first implemented by White 28 ("state prediction") is 
routinely used to predict the outcome of large sparse matrix diagonalizations, which 
no longer occur during time evolution. Here, it merely serves as a basis transformation. 
We will outline the calculation for shifting the active bond by one site to the left. 

Starting from 

1^) = E E E E ^™f_,.,.,..™f,j<-.)i'-^)i'^'+i)i<..)^ (57) 

one inserts the identity X^m^ \''^mi) {'^mi \ obtained from the Schmidt decomposition 
(i.e. density matrix diagonalization) to obtain 

IV') = E EE'^^f-.-^'-f K-i)i^')iO' (58) 
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where 

After inserting in a second step the identity X^mf ^ai^i \'^mi-.2'^i-i) {'^fni-2^i-^\^ ^"^^ 
ends up with the wave function in the shifted bond representation: 

where 

6. Adaptive time-dependent DMRG 

The adaptive time-dependent DMRG algorithm which incorporates the TEBD 
simulation algorithm in the DMRG framework is now set up as follows (details on 
the finite-system algorithm can be found in Ref. 

0. Set up a conventional finite-system DMRG algorithm with state prediction using 
the Hamiltonian at time t = 0, H{0), to determine the ground state of some 
system of length L using effective block Hilbert spaces of dimension M. At the end 
of this stage of the algorithm, we have for blocks of all sizes I reduced orthonormal 
bases spanned by states \mi), which are characterized by good quantum numbers. 
Also, we have all reduced basis transformations, corresponding to the matrices 
A. 

1. For each Trotter time step, use the finite-system DMRG algorithm to run one 
sweep with the following modifications: 

i) For each even bond apply the local time evolution U at the bond formed by 
the free sites to {tp). This is a very fast operation compared to determining 
the ground state, which is usually done instead in the finite-system algorithm. 

ii) As always, perform a DMRG truncation at each step of the finite-system 
algorithm, hence 0{L) times. 

(iii) Use White's prediction method to shift the free sites by one. 

2. In the reverse direction, apply step (i) to all odd bonds. 

3. As in standard finite-system DMRG evaluate operators when desired at the end 
of some time steps. Note that there is no need to generate these operators at all 
those time steps where no operator evaluation is desired, which will, due to the 
small Trotter time step, be the overwhelming majority of steps. 

The calculation time of adaptive time-dependent DMRG scales linearly in L, 
as opposed to the static time-dependent DMRG which does not depend on L. The 
diagonalization of the density matrices (Schmidt decomposition) scales as N^^^^M^; 
the preparation of the local time evolution operator as N^^^.^, but this may have to be 
done only rarely e.g. for discontinuous changes of interaction parameters. Carrying out 
the local time evolution scales as N^^^^^A'P; the basis transformation scales as N^^^^M^. 
As M 3> A'^site typically, the algorithm is of order 0{LN^^^^AI^) at each time step. 
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Figure 8. Time evolution of tlie real part of nearest-neighbor correlations in a 
Bose-Hubbard model with instantaneous change of interaction strength using the 
adaptive time-dependent DMRG. The different curves for different M are shifted 
(parameters as in section 1^. 

7. Case study: time-dependent Bose-Hubbard model 

In this section we present some results of calculations on the Bose-Hubbard 
Hamiltonian introduced in section |21 which have been carried out, using modest 
computational resources and an unoptimized code (this concerns in particular the 
operations on complex matrices and vectors). In the following, Trotter time steps 
down to 5t = 5 X 10~^ in units of h/ J were chosen. It is also important to note that 
in contrast to the DMRG calculations shown earlier for conventional time-dependent 
DMRG up to iVsitc = 14 states per site were used as a local site basis for all calculations 
in this Section. 

Comparing the results of the adaptive time-dependent DMRG for the Bose- 
Hubbard model with the parameters chosen as in section |5] with the static time- 
dependent DMRG wc find that the convergence in M is much faster, for the nearest 
neighbor correlations it sets in at about M = 40 (Fig. |SJ) compared to M — 100 for 
the static method (Fig. |3J). 

This faster convergence in M enables us to study larger systems than with 
static time-dependent DMRG (Fig. O. In the L — 2>2 system considered here, we 
encountered severe convergence problems using static time-dependent DMRG. By 
contrast, in the new approach convergence sets in for M well below 100, which is 
easily accessible numerically. Let us remark that the number M of states which have 
to be kept does certainly vary with the exact parameters chosen, depending if the state 
can be approximated well by matrix product states of a low dimension. At least in the 
case studied here, we found that this dependency is quite weak. We expect (also from 
studying the time evolution of density matrix spectra) that the model dependence of 
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Figure 9. Time evolution of the real part of nearest-neiglibor correlations in 
a Bose-Hubard model with instantaneous change of interaction strength using 
the adaptive time-dependent DMRG but for a larger system L = 32 with Af = 32 
bosons. The different curves for different M are shifted, comparing M = 30, 50, 70 
to M = 80 respectively. 



M is roughly similar as in the static case. 

Similar observations are made both for local occupancy (a simpler quantity than 
nearest-neighbor correlations) and longer-ranged correlations (where we expect less 
precision). Moving back to the parameter set of section |21 we find as expected that 
the result for the local occupancy (Fig. llOfl is converged for the same M leading 
to convergence in the nearest- neighbor correlations. In contrast, if we consider the 
correlation (6^6) between sites further apart from each other the numerical results 
converge more slowly under an increase of AI than the almost local quantities. This 
can be seen in Fig. where the results for M = 40 and M = 50 still differ a bit for 
times larger than t « 2h/J. 

The controlling feature of DMRG is the density matrix formed at each DMRG 
step - the decay of the density-matrix eigenvalue spectrum and the truncated weight 
(i.e. the sum of all eigenvalues whose eigenvectors are not retained in the block bases) 
control its precision. In the discarded weight for the Bose-Hubbard model of section 
12 shown in Fig. ^1 we can observe that the discarded weight shrinks drastically, 
going from M — 20 to M — 50. This supports the idea that the system shows a fast 
convergence in M. Even more importantly, the discarded weight grows in time, as 
the state that was originally a ground state at t < decays into a superposition of 
many eigenstates of the system at t > 0. However, in particular for larger M, it stays 
remarkably small throughout the simulation, indicating that adaptive time-dependent 
DMRG tracks the time-evolving state with high precision. Moving to the detailed 
spectrum of the density matrix (shown in Fig. 1131 for the left density matrix when 
the chain is symmetrically decomposed into S and E), the corresponding distribution 
of the eigenvalues can be seen to be approximately exponential. In agreement with 
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Figure 10. Time evolution of tlie occupancy of tlie second site. Parameters as 
used in section|2]{i = 8, N = 8). The different curves for different M are shifted. 
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Figure 11. Time evolution of the real part of the correlation between site 2 and 
7. Parameters as used in section I^J with N = 8 particles. The different curves 
for different M are shifted. Note that the plot starts at 4 = 1 (parameters were 
changed at i = 0). 
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Figure 13. Eigenvalue spectrum of the left reduced density matrix at different 
times for a symmetric S/E decomposition. Parameters chosen as in section 1^ 
M = 50 states retained. 



the increasing truncation error, one also observes that the decay becomes less steep as 
time grows. Yet, we still find a comparatively fast decay of the eigenvalue spectrum at 
all times, necessary to ensure the applicability of TEBD and adaptive time-dependent 
DMRG respectively. 
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Figure 14. Convergence in the Trotter time of the real part of the nearest- 
neighbor correlations between site 2 and 3 in a Bose-Hubbard model with 
instantaneous change with the parameters chosen as in section |2] at a fixed time. 

Note for all results shown that the unusually large number of states per site 
{Nsite = 14) which would not occur in Hubbard or Heisenberg models could there be 
translated directly into longer chains or larger state spaces (larger M) for the same 
computational effort, given that the algorithm is 0{LN^^^.^M^). In that sense, we 
have been discussing an algorithmically hard case, but in fermionic models DMRG 
experience tells us that M has to be taken much larger in fermionic systems. For the 
fermionic Hubbard model, with A^site = 4, more than A/ = 300 is feasible with the 
unoptimized code, and much higher M values would be possible if optimizations were 
carried out. This should be enough to have quantitatively reliable time-evolutions for 
fermionic chains, while of course not reaching the extreme precision one is used to 
in DMRG for the static case. As the algorithmic cost is dominated by (A^sitc-^^)"^, 
the product NgitcM is an important quantity to look at: while current TEBD 
implementations range at 100 or less, adaptive time-dependent DMRG using good 
quantum numbers runs at the order of 1000 (and more). 

Let us conclude this section by pointing out that at least one improvement can be 
incorporated almost trivially into this most simple version of adaptive time-dependent 
DMRG. Since we have used a first-order Trotter decomposition, we expect that for 
fixed M results of measurements at a fixed time converge linearly with respect to the 
time step dt chosen, as the error per time step scales as St'^, but the number of time 
steps needed to reach the fixed time grows as St~^. In other words, the Trotter error 
is inversely proportional to the calculation time spent. This can indeed be observed 
in results such as presented in Fig. 

It is very easily and at hardly any algorithmic cost that a second order Trotter 
decomposition can be implemented, leading to errors of order St'^. The second order 
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Trotter decomposition reads pPj 

where we have grouped all local Hamiltonians on odd and even bonds into Hodd and 
Heven respectively. At first sight this seems to indicate that at each Trotter time 
step three (instead of two) moves ("zips") through the chain have to be carried out. 
However, in many applications at the end of most time steps, the Hamiltonian does 
not change, such that for almost all time steps, we can contract the second e"'^'"*'*''*/^ 
from the previous and the first Q~^HaddSt/2 fj-Qj-^ i]^^ current time step to a standard 
g-iHoddSt ^jjjig step. Hence, we incur almost no algorithmic cost. This is also standard 
practice in Quantum Monte Carlo 01]; following QMC, second order Trotter evolution 
is set up as follows: 

(i) Start with a half-time step e-i^odd<5t/2_ 

(ii) Carry out successive time steps e^'^"""""^* and e^'^"'*''''*. 

(iii) At measuring times, measure expectation values after a e~^^°'''^^* time step, and 
again after a time step 6"'^""="**, and form the average of the two values as the 
outcome of the measurement. 

(iv) At times when the Hamiltonian changes, do not contract two half-time steps into 
one time step. 

In this way, additional algorithmic cost is only incurred at the (in many applications 
rare) times when the Hamiltonian changes while strongly reducing the Trotter 
decomposition error. Even more precise, but now at an algorithmic cost of factor 
5 over the first or second-order decompositions, would be the usage of fourth-order 
Trotter decompositions (leading to 15 zips through the chain per time step, of which 
5, however, can typically be eliminated) 03 ESI ■ 



8. Conclusion 



The TEBD algorithm for the simulation of slightly entangled quantum systems, such 
as quantum spin chains and other one-dimensional quantum systems, was originally 
developed in order to establish a link between the computational potential of quantum 
systems and their degree of entanglement, and serves therefore as a good example of 
how concepts and tools from quantum information science can influence other areas 
of research, in this case quantum many-body physics. 

While exporting ideas from one field of knowledge to another may appear as an 
exciting and often fruitful enterprise, differences in language and background between 
researchers in so far separated fields can also often become a serious obstacle to the 
proper propagation and full assimilation of such ideas. In this paper we have translated 
the TEBD algorithm into the language of matrix product states. This language is 
a natural choice to express the DMRG algorithm - which, for over a decade, has 
dominated the simulation of one-dimensional quantum many-body systems. In this 
way, we have made the TEBD algorithm fully accessible to the DMRG community. 
On the other hand, this translation has made evident that the TEBD and the DMRG 
algorithms have a number of common features, a fact that can be exploited. 

We have demonstrated that a very straightforward modification of existing finite- 
system DMRG codes to incorporate the TEBD leads to a new adaptive time-dependent 
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DMRG algorithm. Even without attempting to reach the computationaUy most 
efficient incorporation of the TEBD algorithm into DMRG implementations, the 
resulting code seems to perform systematically better than static time-dependent 
DMRG codes at very reasonable numerical cost, converging for much smaller state 
spaces, as they change in time to track the actual state of the system. On the other 
hand, while it presents no new conceptual idea, the new code is also significantly 
more efficient than existing embodiments of the TEBD, for instance thanks to the 
way DMRG handles good quantum numbers. While we have considered bosons as an 
example, as in standard DMRG fermionic and spin systems present no additional 
difficulties. Various simple further improvements are feasible, and we think that 
adaptive time-dependent DMRG can be applied not only to problems with explicitly 
time-dependent Hamiltonians, but also to problems where the quantum state changes 
strongly in time, such as in systems where the initial quantum state is far from 
equilibrium. The method should thus also be of great use in the fields of transport 
and driven dissipative quantum systems. 
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